function [ energy ] = trapEnergy( x,y,z )
% computes the trap configuration energy

energy = 0;
for i = 1:length(x)
    for j = i+1:length(x)
        s = sqrt((x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2);
        energy = energy + 1/s - 1/2*(log(s) + log(2+s));
    end
end

end

